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r-| ' We study a system of A identical interacting bosons trapped by an external field by 
solving ab initio the many-body Schrodinger equation. A complete solution by using, for 
d \ example, the traditional hyperspherical harmonics (HH) basis develops serious practical 
O^' problems due to the large degeneracy of HH basis. Symmetrization of the wave func- 
- ■— j ■ tion, calculation of the matrix elements, etc., become an immensely formidable task as 
^ , A increases. Instead of the HH basis, here we use a new basis, called "potential harmon- 
ics" (PH) basis, which is a subset of HH basis. We assume that the contribution to the 
orbital and grand orbital [in 3(^4 — l)-dimensional space of the reduced motion] quan- 
tum numbers comes only from the interacting pair. This implies inclusion of two-body 
correlations only and disregard of all higher-body correlations. Such an assumption is 
ideally suited for the Bose-Einstein condensate (BEC), which is required - for experimen- 
tal realization of BEC - to be extremely dilute. Hence three and higher-body collisions 
are almost totally absent. Unlike the (3 A — 4) hyperspherical variables in HH basis, the 
PH basis involves only three active variables, corresponding to three quantum numbers 
- the orbital Z, azimuthal m, and the grand orbital 2K + I quantum numbers for any 
arbitrary A. It drastically reduces the number of coupled equations and calculation of 



the potential matrix becomes tremendously simplified, as it involves integrals over only 
three variables for any A. One can easily incorporate realistic atom-atom interactions in 
a straight forward manner. Wc study the ground and excited state properties of the con- 
densate for both attractive and repulsive interactions for various particle number. The 
ground state properties are compared with those calculated from the Gross-Pitaevskii 
(GP) equation. We notice that our many-body results converge towards the mean field 
results as the particle number increases. 

FAGS number(s): 03.65.Ge, 03.75.Hh, 03.75.Nt, 31.15.Ja 

Key words: Bose Einstein Condensation, Hyperspherical harmonics method. Potential 
harmonics. 
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I. Introduction 



Although the phenomenon of Bose Einstein Condensation (BEC) was 
known for a long time [1-3], its experimental observation in trapped and supercooled 
(down to nano Kelvin temperatures) alkali atoms in 1995 [4-6] renewed a great deal of 
interest - both experimental and theoretical - in the phenomemon. The importance of 
this topic is clearly demonstrated by the fact that two independent Nobel Prizes were 
awarded on BEC related works in quick succession in the recent past. The density 
of magneto-optically trapped atomic gas undergoing BEC is extremely low ( to avoid 
recombination of atoms through three and higher body collissions) and the number of 
trapped atoms is typically of the order of a few hundred to a few million. This is 
extremely small compared to the Avogadro number. For such a small number of atoms 
an exact ab initio solution would have been ideally desirable. But an interacting system 
oi A = {N + 1) particles has 3N relative degrees of freedom and an ab initio solution of 
the corresponding Schrodinger equation is practically impossible for A > 3. Hence the 
usual theoretical tools that have been used so far are the mean field models [7-10] and the 
Thomas-Fermi [8] approximation. The dilute atomic gas undergoes BEC below a critical 
temperature ( typically 10~^ degree K) when most of the atoms (bosons) go to the 
single particle ground state. Then the de Broglie wavelength associated with the atomic 
motion is much larger than the interaction length scale. Hence the resulting many body 
system emerges as essentially a single quantum system where all the atoms behave in a 
coherent manner [8,11]. At zero temperature, the effect of the excited states are absent 
and the condensate is described by a single equation involving the condensate wave 
function [8]. However this simple picture is no more true at a finite temperature due to 
the existence of interparticle interactions. The usual procedure is to start with the mean 
field approximation fike the Hartree-Fock (HF) theory for the many body system [7-10]. 
This is an independent particle approach where each individual atom is assumed to move 
in a single particle orbit. These orbits are determined self consistently by allowing an 
atom in one orbital to be infiuenced by other atoms in other orbitals through two-body 
interaction. Assuming a contact interaction for the two-body potential, viz., V{f—r') — 
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g5(f—r'), the many body equation reduces to the famous Gross-Pitaevskii (GP) equation 
[8]. At zero temperature, the effect of excited states are neglected and the condensate is 
described by the time independent GP equation 



^ v'+V;xt(r) + <70'(f) 
2m 



0(f) = /i0(f) (1) 



where n(r) = is the condensate density and /i is the chemical potential. For a first 
approach the contact interaction is justified since in the cold and dilute gas only binary 
coUissions at low energies arc relevant. These arc characterized by the s-wavc scattering 
length (osc), which is independent of the details of two-body potentials. The strength 
constant g of the contact interaction is related to the scattering length through [8] 

9 = 2 

m 

The GP equation has been used extensively to study the BEG [8,11]. Although most 
of the static, dynamic and thermodynamic properties are fairly well reproduced by the 
GP equation [8], the wave function does not include any correlation. Furthermore the 
assumption of a contact 5-interaction is too simple and does not represent the realis- 
tic situation. It has already been shown that the Dirac (5-function is not suitable as a 
replacement of the actual two-body interaction in exact theories in more than one dimen- 
sion [12]. This is because the Hamiltonian then becomes unbound from below and the 
ground state energy diverges for an attractive zero range potential. Solutions are usually 
obtained in the metastable region (although such solutions are not rigorously correct for 
an attractive (5-function potential) and the condensate becomes unstable for N larger 
than a critical number, due to disappearance of the local minimum. This was shown by 
Bohn et al in a hyperspherical calculation keeping the lowest (most dominant) harmonic 
[13]. A third disadvantage is the non-lineraity of the GP equation, so that standard 
quantum mechanics is not applicable without concessional approximation. Thus one has 
to go beyond the mean field approximation and simple contact interactions. 

Because of the limitations of the mean field theory and GP equation it 
is desirable to solve the many body linear Schrodinger equation directly. The Schrodinger 
equation for a system oi A — {N identical bosons, each of mass m, confined by an 



external field T^^ap (acting on each individual boson) and interacting through a mutual 
two body interaction V is 

-^2 A A A 

2m 



i=l i=l i<j=2 



= E^{x) (3) 



where x refers to the set of particle coordinates {xi,X2, xa} of A bosons. The cen- 
ter of mass (CM) motion can be eleminated resulting in a Schrodinger equation in 3A^ 
variables. A standard practice is the use of hyperspherical harmonics expansion (HHE) 
method, in which the wave function is expanded in the complete set of hyperspherical 
harmonics (HH) spanning the (3A^ — l)-dimensional hyperangular space [14]. Projection 
on a particular HH leads to a system of coupled differential equations (CDE). However 
there are several very serious difficulties associated with the solution of a fairly large num- 
ber of particles. Firstly the expansion basis of HH should be properly symmetrized and 
appropriate conserved quantum numbers properly taken care of. Secondly calculation of 
matrix elements of all the pairwise two-body potentials is an extremely formidable task. 
Finally, due to very large degeneracy of the HH basis for a large number of particles, the 
number of CDE and the dimemsion of the potential matrix is too large to be handled 
by any computer [14]. On top of all these, the convergence rate of the HH expansion, 
especially for long-range interactions, is slow [15]. For these reasons the HHE method 
has been used fully for the three body system only [15-18]. On the other hand, as we 
discussed earlier, the condensate can be treated broadly as a "single lump of quantum 
stuff' , since all the individual atoms in the condensate lie within one single de Broglie 
wavelength [8]. Thus it is reasonable to assume that the basic properties of the conden- 
sate in the lowest approximation, is described by a single collective coordinate. This led 
Bohn et. al. [13] to go for the K-harmonic approximation, in which the HH expansion 
is restricted effectively to the first term only ( which is independent of the hyperangles) . 
Such a drastic approximation may be justified for a contact interaction only. Even in 
this case, for an attractive 5-function interaction, there are no rigorously stable solutions. 
Since the wave function becomes independent of the hyperangles and the hyperradius 
is invariant under any permutation of the particles, the wave function becomes totally 
symmetric, as required. The calculation of the potential matrix also simplifies immensely 



and the CDE reduces to a single differential equation [13]. The hyperradius emerges as 
the sought for collective coordinate. In spite of the great simplifications, there are serious 
criticisms of this approach : (1) The method cannot be applied to any realistic two-body 
interaction. (2) Even for a contact interaction, the method is not satisfactory for attrac- 
tive (5-f unction interaction, for which no rigorous solution extsits. (3) Only one collective 
variable is involved. Hence it can only describe the gross features of the condensate, 
without any finer details. Thus a more rigorous treatment is necessary. But as already 
mentioned a completely rigorous, essentially exact solution of the Schrodinger equation is 
possible for the three body system only. That has been done to get an idea of the initial 
trend as the particle number increases from three by Esry and Greene [12]. However that 
is far from the real situation in a condensate. 

An alternative approach of exact numerical diagonalization of the many 
body Hamiltonian was adopted by Haugset and Haugerud [19] for a small number (< 30) 
of interacting (via contact interaction) bosons confined by a harmonic trap. However, this 
was restricted to one and two dimensions only. Moreover the process is extremely time 
consuming even for two dimensional condensates, with a nagging question of convergence 
of the chosen harmonic oscillator basis expansion. The rate of convergence is expected 
to be slower for a realistic two-body interaction and in three dimensional condensates. 
Although analytic expressions for the matrix elements are greatly simplified for a delta 
function interaction, all the problems associated with a contact interaction discussed 
above remain for the two dimensional condensate. However, there is no problem with 
the one dimensional condensate, as one dimensional delta function is not pathological. 

From the above discussion it is clear that an exact treatment of the 
many body system in three dimensions is not possible beyond the three body system. 
On the other hand the single quantum nature of the entire condensate suggests that out 
of the thousands to millions of degrees of freedom of the individual particles only a few 
arc physically relevant. This is due to the fact that the condensate is possible only at 
extremely low temperatures ( low energy of the individual particles) and extremely low 
densities. Under these conditions only two body collisions are relevant. Three and higher 
body coUsions are extremely rare and correlations beyond two body correlations in the 
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condensate wave function are completely negligible upto a very high degree of precision. 
Indeed in an experimental situation this is ensured by keeping the density extremely 
low, so that there are no recombination via three and higher body collisions [8]. The 
mean field approach ignores all correlations including two-body correlations. Importance 
of two-body correlations in BEC has been emphasized by several authors [20,21]. Thus 
physically relevant quantities are contributed by two-body collisions, while the rest of the 
particles in the condensate do not partate in any motion other than a collective one and 
are simply inert spectators. The emerging picture then suggests that most of the degrees 
of freedom of these spectators can be frozen, while a single pair interacts. This reduces 
the physically important degrees of freedom of the condensate to just four - a global 
length scale (hyperradius) of the entire condensate, and the three degrees of freedom of 
the relative vector r^j = Xi — Xj of the interacting pair. However one has to concede that 
any pair out of the A = {N + 1) atoms in the condensate can interact. These are also 
consistent with the intuitive "single quantum stuff' concept of the condensate. 

Among the various possible theoretical approaches to handle the many 
body system, the HHE method appears to be the most lucrative one, as it readily pro- 
vides the hyperradius as the most important collective variable. A theoretical formahsm, 
arising out of the HHE method, was adopted by Fabrc de la Ripelle [22] in 1986. Al- 
though the primary concern there was an application to the nuclear systems consisting of 
fermions, it was noted that the formalism is applicable to a system of identical bosons also 
[23]. To incorporate the importance of the interacting pair and two-body correlations, 
he introduced the potential harmonics (PH) expansion basis [23] , rather than the general 
HH basis, thereby reducing the expansion basis to a great extent. Potential harmonics 
is a subset of HH, where all correlations higher than two-body ones are disregarded. In 
PH, the contribution to the total orbital angular momentum as also the grand orbital 
quantum number comes only from the interacting pair. Here all the (^ — 2) spectators are 
assumed to be described by the HH of the lowest (zero) order. We adopt this procedure 
since this approximation is quite justified in our situation due to the diluteness of BEC, 
where two-body correlation is the most important and all higher-body correlations can 
be safely ignored. Using Faddeev hke decomposition of the total wave function, and then 
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expanding each such component in an appropriate set of PH, the number of CDE can be 
reduced drastically. Since the PH involves only four active degrees of freedom, calcula- 
tion of potential matrix elements is simplified tremendously as compared to that in HH 
basis. Use of realistic two-body interactions and calculation of their matrix elements are 
quite straight forward. Requiring the Faddeev component for the (ij) interacting pair to 
be symmetric under (i_7)-pair exchange, the total wave function becomes automatically 
totally symmetric. Thus the symmetrization of the wave function is also handled prop- 
erly. 

Thus a truely many body equation is reduced to a tractable mathe- 
matical form. The assumptions leading to this are especially appropriate for the EEC. 
Hence we adopt the PH basis as our starting point. This is theoretically applicable to 
a system containing any number of particles, but we will see in Sec. HI, that numerical 
difficulties arise as the number of particles increases beyond a certain number. In this 
communication we report some of the basic properties of the condensate for various par- 
ticle numbers and compare them with previous calculations. 

Sorensen et al [20,21] have followed a method which is similar in spirit 
to the present work, although it differs in details. They expand the wave function in the 
adiabatic subset ^n{p,^) of the full {N — l)-body Hamiltonian (in CM frame). Later 
this is decomposed in Faddeev like components (pij. This leads to an integro-differential 
equation (IDE) for {—(t>ij-i which is the same for all ij-pairs due to boson symme- 
try) involving five dimensional integrals and the full (3A^ — 4)-dimensional hyperangular 
differential operator A^. All (3iV — 5) angle derivatives other than a = ayi (where 
= \/2psinaij is the relative separation of the (ij)-pair and p is the hyperradius of 
the full system) , are disregarded, leaving only one angle variable. Assumption of a very 
short ranged two-body potential reduces the five dimensional integrals to two dimen- 
sional ones. In this limit simple expressions are obtained for the integrals in IDE. On the 
other hand, we write the complete 3A^-dimensional Schodinger equation of the relative 
motion of a (A^ + 1) boson system in terms of Faddeev like components $(r7j, r), subject 
to the approximation that $(r7j , r) corresponds to zero eigenvalue of the hyper angular 
momentum operator (see later) for the (A^ — 1) remaining relative vectors of the specta- 
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tors, while (ij)-pair interacts. These are then expanded in the potential harmonics (PH) 
basis. The assumptions in our method are clearly justified in terms of the physics of the 
chosen system, which have been stated earlier. While the use of PH basis in nuclei (as 
originally used by Fabre in [22,23]) is questionable due to high spatial density of nucleons 
in a nucleus, its application in BEC is ideally suited (the number of atoms in the con- 
densate is < 10^ in a space of macroscopic linear dimensions of order 10~^cm, which is 
immensely smaller than the Avogadro number). As a consequence, the total orbital (/) 
and grand orbital (K) angular momenta of the system are contributed by the interacting 
pair alone. Apart from this well justified fundamental approximation, we need no other 
approximation. Although for the first calculation, we have restricted ourselves to Z = 
and a central two-body interaction, both these can be relaxed resulting in a somewhat 
more complicated equation. Finally the system of coupled differential equations in one 
variable (hypcrradius, r) can be solved numerically, without additional approximation (as 
done in ref. [15] and compared with adiabatic approximation (AA) in ref. [28]) using, 
e.g. renormalized Numerov method. Once again, as a preliminary calculation, we use 
AA to solve the CDE. Our use of AA in solving the CDE is not an indispensable one; 
it is done only to reduce the numerical calculation. But in the approach of Sorensen et 
al, adiabatic subset is the starting point to separate the hypcrangular and hyperradial 
motions. Furthermore our method can handle any two-body potential (central or not, 
short ranged or not) ; for non-central potential, calculation of matrix elements will involve 
integrals over two polar angles in addition. The approach of ref. [20,21] requries a very 
short ranged, central potential to reduce the equation to a manageable form. The present 
method has no such restriction. 

The paper is organised as follows. In Sec. II, we present our choice 
of Jacobi coordinates and express the kinetic energy in the chosen set of hyperspherical 
variables. In the same section, we introduce the concept of potential harmonics basis and 
obtain the set of coupled differential equations resulting from the many-body Schrodinger 
equation. The numerical method for solving the CDE and results of our calculation are 
presented in Sec. III. There we compare our results for different numbers of particles 
with those of earher calculations. Finally in Sec. IV we draw our conclusions. Some of 
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the detailed expressions have been given in the Appendix. 



II. Theory 

A. Choice of Jacobi coordinates 

We consider a system oi A — {N + 1) identical bosons, each of mass m 
and confined magnetically in a trap which is approximated by a spherically symmetric 
harmonic oscillator potential with frequency uj. We assume that the atomic cloud is at 
zero temperature. The full many body Hamiltonian is given by 



^2 7V+1 N+1 ^ N+1 

^^l^ i=l i=l ^ ij>i 



= E'^{x) (4) 



where x refers to the set of particle coordinates {xi, X2, x^+i} of {N + 1) bosons and 
E' is the total energy. We decompose the total wave function ^{x) as the sum of pairwise 
partial waves 

N+1 
ij>i 

The Schrodinger equation for ipy can be written as 

N+l 

(T + K.ap - E')AA^) = -y (n,) M^) (6) 

kl>k 

where T is the total kinetic energy operator, V/^^^ is the confining potential, V^^^p — 
Y^^Ji^ \muj'^xl and V{rij) is the pairwise local central two-body interaction between i*'* 
and particles, flj = Xi — Xj. Applying the operator J2ijti both sides of eq.(6), 
and using eq.(5), we get back eq.(4). Now instead of (A^-|- 1) particle coordinates Xi, the 

— * 

system can alternatively be described by the center of mass coordinate R 

I ^+1 



R 



N + 

and N Jacobi coordinates defined as 



Y E ^1 (7) 
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The chosen normahzation of facihtates writing the Laplace operator in the form 

N 

Til. -\- 

2 i^, 2A 



1 1 



(9) 



i=l i=l 

Then the relative motion (after removal of center of mass motion from eq.(4)) is described 
by [14,23] 



2 AT 



^ 1=1 



V'(Ci,-,Civ)=o , 



(10) 



- N 


1/2 


2 


1/2 


Ecf 

j=i . 









where Vtrap — Z^j=i f C ^^id ^nt is the sum of all pairwise interactions, Vint — 
^(^y ) expressed in the relative coordinates. Here E is the energy of the relative 
motion, i.e., E' minus energy of CM motion. The hyperradius r is defined as [22] 



(11) 



which is invariant under permutations of the particle indices as also three dimensional 
rotations. The hyperspherical coordinates are constituted by the hyperradius r and 
remaining (3A^ — 1) hyperangles, denoted coUectivelty by Qat in -D = 3A^ dimensional 
space. Note that the choice of Jacobi coordinates eq.(8), is not unique, since the labelling 
of the particle indices and consequently that of the Jacobi coordinates are arbitrary. We 
choose a particular set by specifying the relative separation of the interacting pair, r7j as 
Civ and {'&,(p) are the two spherical polar coordinates associated with fij. The relative 
length is defined in terms of through r^^ = r coscf). For the rest of {N — 1) Jacobi 
coordinates, we define the hyperradius pij in the 3(A'' — 1) dimensional space by 



Af-l 



Pij 



:E4T/^ 



(12) 



k=l 



Hi] I ij 



which is related with Qjsi = r\j by 

^, Pij — r sirKp 
Then our hyperspherical coordinates become 

(r, Qat) = (r, (p, (fi, D.N-i) 



(13) 



(14) 
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Here JIat-i involves 2(A'" — 1) spherical polar angles associated with each of (N — l) Jacobi 
vectors {Ci, C2, Ca^-i} and {N — 2) angles (expressing relative lengths) , i.e., a total of 
(3A^ — 4) variables. In this choice of hyperspherical coordinates, the Laplace operator 
takes the form [22] 

V^-Evl = ^ + ^^ + ^, A = N + 1 (15) 

LP'iyt]^) is the grand orbital operator in ?>N dimensional space which is obtained from a 
recurrence formula [22] and has the form 

L^(Q^) = 4(1 - .^)|^ + 6[2 - N{1 + z)]^ + 2^ + 2^lf^ (16) 
oz-^ oz 1 + z 1 — z 

where z — cos2(l), cuij reprsents the two polar angles (p) associated with r7j and 

L'^{Qn-i) is the grand orbital operator in 3(A^ — 1) dimensional space. 

B. Potential basis and potential multipoles 

To exapand a function V{rij) in hyperspherical harmonics (HH) we 
use the above definition of Jacobi coordinates. It is easy to sec that HH basis which is 
complete for the expansion of V{rij) does not contain any function of the coordinate Q 
with i < N and is given by [23] 

r'2K+li^^J) = '^''^P2'K+li<P)yo{D - 3) (17) 

where ^""^P^k+i is a function involving the Jacobi polynomial and is needed in the general 
expression of the hyperspherical harmonics (see Appendix) of grand orbital 2K + I and 
orbital angular momemtum /. The quantity yo{D — 3) is the HH of order zero (i.e. 
grand orbital quantum number is zero) in 3(A'^ — 1) dimensional space, yo{D — 3) = 
( J^^lfo'-a)/^^ ) ^ ■ '^^i^ basis set which is a subset constituted by HH of order {2K + 1) 
are called "potential harmonics" (PH). These are the eigenfunctions of L^(f2Ar), when 
the eigenvalue of L^(r2jv-i) is and satisfy the eigenvalue equation : 

L^{nN)+jC{C + D-2)]vl^K+i{%)^0, C^2K + l ■ (18) 

The relation L^(Oiv-i)'0ij(^) = implies that we are considering only those states which 
are invariant under all generalized rotations in 3(A'" — 1) dimensional space. Natuarally 
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the contribution to the grand orbital quantum number comes only from the interacting 
pair. This corresponds effectively to two-body correlations only in the wave function. 
Due to diluteness of atomic BEC, the effect of higher body correlations can be ignored 
as the probability for three or more particles to come close at the same time is extremely 
small. This reduces the number of quantum numbers in the new basis (all the quantum 
numbers specifying the eigenfunctions of L^(Qjv_i) are zero). It contains only three 
quantum numbers; orbital I, azimutal m and grand orbital 2K + 1 for any N, instead of 
{3N — 1) quantum numbers corresponding to (3A^ — 1) hyperspherical variables in the 
general HH basis. The normalization condition is given by 

j 'P2K+l{^ij)'P2K'+l'{^ij)d^ij = ^KK'^ll'^mm' ' (19) 

Then the PH expansion of the potential is 

V{nj) = E ^r(^,j)^2'K+^(%)4^''^(r) (20) 

K,l,m 

Aj"(i, j) is an operator which is independent of r^j, but may act on other variables like 
spin variables. The quantity vj;P'^\r) are the "potential multipoles" and for a central 
potential, it is given by [23] 

= |3^o(^-3)|-Vo"/'(^)^2^+^(0)l^(rW)(sm0)^-4(W)2# , 

where the functions '^^^ P2k+i{^) defined in the Appendix. Starting from the multi- 
poles calculated either for the D = 5orD = 6( depending wheather D is odd or even) 
and using simple recurrence formulae potential multipoles for any D can be calculated 
[23]. 

C. Coupled differential equations 

Splitting eq.(lO) in the manner of eq.(6) for the (?j) -interacting pair and using eqs.(14)- 
(16), subject to the restriction that the eigenvalue of L^(f2Ar_i) is zero, wc see that the 
{ij) Faddeev component will be a function of r^j and r only and satisfies [23] 

(T + Vtrap - ^)$(r^-, r) = -V{n,) ^ $(rw, r) , (22) 

k,l>k 
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where ^{fij,r) differs from the general solution ipij by the fact that it corresponds to 
eigenvalue zero of the operator L'^^Q^^i). Next expand the wave function $(ri^ , r) in the 
complete set of potential harmonics (when / is a good quantum number) as 

mj,r)=r-'^j:'P2K'+i{^rjWAr) (23) 

K' 

Substitution of eq.(23) in eq.(22) and projection on the same basis, leads to the set of 
coupled differential equations [23] 



(f Ck{Ck + 1) 



where 



+ ^ V ' +VtraAr)-E 



u^Kir) + E fl'iVKArWAr) = (24) 

K' 



fxi — T,k,l>k < 'P2K+l(^ij)\'P2K+l(^kl) > 

The potential matrix is given by 

VKK'ir) = I V'^\i{n,^)V{n^)P^,^,{n,^)dnr, (26) 

So instead of (3A^ — 1) angle variables in HHE method, in potential harmonics expansion 
method (PHEM) the integral invloves only 3 angle variables. It greatly simplifies the 
calculation of the matrix element for any N. 

The quantity /|; of eqs. (24) and (25) is given by 



/i, = 1 + [2{A - 2)(-l)'P«^(-l) + 2)^(-^ ^^ Pf{-l)5iAlPf{l) (27) 

where a — {2>A — 8)/2 and (5 — I -\- ^ and P^^{x) is the Jacobi polynomial. Multiplying 
eq. (24) by appropriate constant factors, it can be put in a symmetric form: 



(28) 



+ Y.K'VKK'{r)UK'i{r)^Q 
where C — I -\- {?>A — 6)/2, the symmetrized potential matrix V kk' has the form 

VKK'{r) = fKiVKK'{r)fK'i{hfhf,)-"^ (29) 

and 

UKiir) ^ fKiihff^u'^ir) ■ (30) 
14 



Here hf is the norm of the Jacobi polynomial P^^{x) [24]. The potential matrix element 
is obtained from eq. (26), using eq. (17) and eq. (42) of Appendix, in the form 

VKK'ir) = Pf{z)V (1^/^) P^^{z)w,{z)dz, (31) 

where wi{z) = {1 — z)"-{l + z)'^ is the weight function of the Jacobi polynomials [24]. 
For Gaussian interaction with A = 3, the integral can be obtained analytically [25], from 
where one can directly check the numerical accuracy. 

III. Numerical method and results 

A. Numerical method 

For a chosen number of particles (^4) and a chosen interaction potential 
{V{rij)), we calculate the potential matrix for a fixed value of hyperradius (r) from cqs. 
(29) and (31) using a multi-point Gauss-Jacobi quadrature. For the present calculation 
we select I — and truncate the PH expansion basis of eq.(23) to a maximum K value 
(= Kjnax)- In order to simplify the solution of the set of coupled differential equations, 
eq. (28), we adopt the hyperspherical adiabatic approximation (HAA) [16,26]. In this 
approximation it is assumed that the hyperradial motion is slow compared to the hyper- 
angular motions. Hence the latter can be solved adiabatically for a fixed value of r to get 
an effective potential as a parametric function of r [16]. This is done by diagonalizing the 
potential matrix together with the diagonal hypercentrifugal repulsion and the trapping 
potential for each value of r : 

E MKK'ir) XK'x(r) = u;x(r)xKx(r) (32) 

K'=l 

where 

The lowest eigenvalue gives the "lowest eigen potential ", u!o{r). As we discussed in 
the introduction, the hyperradius behaves as the most important collective coordinate 
and a;o(r) is the potential in which the condensate moves as a "single quantum stuff', 

15 



mr- 



;{C{C + 1) + AK{K + a + + 1)} + Vtrap{r) 



JKK' 



(33) 



except for attractive two-body interactions and A > Acr (see later). Another collective 
coordinate is the hyperangle appearing in the wavefunction through eqs. (23) and 
(17), which describes the deviations of the condensate from hyperspherically symmetric 
distribution. 

In the HAA approach, an approximate solution of eq. (28) is obtained 
by solving a single uncoupled differential equation [16] 

Co(r) = • (34) 

The solution of eq. (34) subject to appropriate boundary conditions on Co('") gives the 
energy E, which is an upper bound for the eigen energy of eq. (28). The partial waves 
of eq.(28) are given in HAA by [16] 

UKi{r) ~ Co(r)xKo(r) (35) 

This approximation is usually called uncoupled adiabatic approximation (UAA) in the 
literature [16,26]; disregarding the third term on the left side of eq.(34) one gets the so 
called extreme adiabatic approximation (EAA). It has been shown that the HAA is in 
very good agreement (having less than 1% error) with the exact solution of the CDE for 
both atomic [27-29] and nuclear [30-31] cases. Since this is adequate for this preliminary 
application of this new method, we adopt the HAA, instead of solving the full set of CDE 
by exact numerical algorithm like the renormalized Numerov method [32] . 

B. Choice of two body interaction potential 

In this report we compare our results with those of the GP equation as 
also with other calculations using a contact (5-interaction. But a (5-function interaction 
is not a physical one since it diverges at rjj = and nothing {e.g. centrifugal repulsion) 
can prevent its overwhelming effect. As a result, the Hamiltonian becomes unbound 
from below for an attractive 5 interaction. This is manifest in the effective potential 
a;o(r), which for a particle number (^4) less than a critical value {Acr) produces a local 
minimum at a finite value of r ( giving rise to a metastable solution), but a;o(r) — > — oo 
as r ^ for any number of particles (see following subsection, as also ref. [12]). Thus 
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+ ^o(r)+ \^ 1—7—1 



m dr"^ dr 



there are no rigorously acceptable and stable solution for any A, since the attractive 
essential singularity at r = will pull the system to r ^ and the corresponding wave 
function will diverge at r = 0. Although the 5-function is particularly convenient for 
analytic calculations, it is desirable to choose an interaction which would either remain 
finite or at worst introduce a removable singularity as rjj — > for attractive cases. Then 
the hyper centrifugal repulsion in eq.(28) (which is non vanishing even for l—O, K — 
and increases rapidly as A increases) will not allow the interacting particles to come too 
close to each other. We thus choose a Gaussian potential of strength Vq and range Tq 



Choosing appropriate values of Vq and ro, the potential can be made either soft or stiff. 
A particular experimental situation at the low temperature limit is characterized by the 
s-wave scattering length (ogc)- For given values of Vq and Tq, one can calculate a^c by 
solving the two-body radial Schrodinger equation for positive energies, in the zero energy 
limit. Alternately, for a suitably chosen value of ro and an experimentally known value of 
asc one can find Vq numerically from the solution of the two-body Schrodinger equation 
in the E 0+ limit. In Fig. 1, we present a plot of calculated function of Vq 

for ro = 0.0855 o.u.. As is well known, Ugc is positive and monotonically continuous for 
Vq > 0. The scattering length becomes negative as becomes negative and continues to 
— oo at a particular negative value of Vq. At this point, a^c has an infinite discontinuity 
and as Vq decreases further, a^c starts from + oo and decreases continuously to —oo at a 
second particular value of Vq. The first, second, branch of the curve (as Vq decreases 
from positive values) correspond respectively to zero, one, two-body bound states. 
For a stable BEC, we choose the first branch of the curve. From Fig. 1, one notices that 
for ro = 0.0855 o.u., the first discontinuity occurs at about Vq = — 184 o.u.. For tq = 
0.005 O.U., this value is much more negative (—8.18963 x 10^ o.u.). In the same figure, 
we also plot the Born approximation for a^c (corresponding to ro = 0.0855 o.u.), given 



V{rij) = VQe < 



(36) 



by [12] 




(37) 
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For a Gaussian interaction this integral can be done analytically and gives 

4? = fv^rl^. (38) 

Prom Fig. 1, it is seen that the Born approximation is good only for small values of |Vo|. 
In this work, we use the exact result and not the Born approximation. For repulsive 
potentials, we choose a conveniently small value of ro and calculate Vq by the exact pro- 
cedure. 

Choosing a smaller value of Tq, Vq increases in magnitude and the po- 
tential becomes stiffer. For very small values of ro, V{rij) simulates a 5-function. For 
attractive interactions, we perform a model calculation with chosen values of ro and Vq. 

C. Results 

With this choice of potential we have solved the CDE eq.(28) for various 
number of particles. We use oscillator units {o.u.) in which energy and length are ex- 
pressed in units of oscillator energy and oscillator length {hcu and ^f^^ respectively, where 
u is the circular frequency of the harmonic confining potential). The matrix element, 
eq. (31), has been calculated by a multi-point Gauss- Jacobi quadrature, the number of 
points being decided by the condition of convergence of a typical matrix element. We 
first verify that our results are independent of the choice of ro, if Vq is appropriately 
calculated using two-body Schrodinger equation, so that asc has the same value ( 100 
Bohr for ^'^Rb, which has a repulsive interaction). In a few representative calculations, 
the ground state energy and low lying excitation spectrum of the condensate containing 
A particles have been found to be stable within numerical errors, for several values of ro 
ranging from 0.1 o.u. to 0.005 o.u. As for example, the ground state energy per particle 
for a condensate containing A — 10 bosons approaches a convergence as ro decreases 
from 0.1 to 0.005. Relative change in the energy per particle from ro — 0.01 o.u. to 0.005 
o.u. is only about 0.012%. As ro decreases, the calculation of the matrix elements as 
also the solution of eq. (34) become extremely CPU time consuming. This is because 
for very small ro, one has to introduce very fine r-mesh intervals ( typically 10~^ o.u.), 
which increases CPU time enormously. To keep the numerical calculations manageable, 
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we choose ro = 0.005 o.u. and Vq — 3.1985 x 10^ o.u. (which corresponds to JILA ^"^Rb 
experiments with a^c = 100 Bohr and trap frequency v = 200 Hz). We next test the 
convergence of our results as Kmax increases by calculating the ground state energy per 
particle of the condensate for Vq — 3.1985 x 10^ o.u. and ro — 0.005 o.u. Our results are 
presented in Table 1. It is seen that the energy per particle converges quite rapidly as 
Kmax increases. For example, for A — 20, the change in energy is less than 0.001% as 
Kmax increases from 2 to 10. Another interesting observation is that the ground state 
energy decreases as Kmax increases, which is consistent with the Rayleigh-Ritz principle. 
Thus it is reassuring that our method is working satisfactorily and is fast converging. 

However a numerical difficulty appears as the particle number {A) and 
Kmax increase. The quantity a increases rapidly with A, (e.g., a — 0.5 for ^4 = 3 and a 
= 71 for A = 50), while (3 remains constant at \ (for / =0). Thus the Jacobi polynomial 
{P^'^{z)) as also its weight function {wi{z)) are highly asymmetric functions in the inter- 
val [-1,1] (see ref. [24]). They have tremendous variation in their values (e.g. 2" to zero) 
as the argument varies from -1 to +1 for large A. Furthermore wi{z) increases from to 
2" within a very small interval close to z — —1, for large a. In addition, P^'^{z) has n 
nodes in the interval [-1,1]. Hence unavoidable numerical error creeps into the numerical 
integration of the potential matrix, using eq. (31). Consequently the calculated energy 
per particle and other physical quantities show irregularity for ^4 > 40, as also for smaller 
A with large Kmax- Therefore we have restricted A to 35. Even for 15 < ^4 < 35, some 
results for large Kmax are not rehable. Hence these have been left out in Table 1. In all 
subsequent calculations, we keep Kmax = 4. We are at present trying to overcome these 
difficulties for large A by improved numerical techniques. 

In Fig. 2, we present a plot of the lowest eigen potential in EAA for A 
— 20, for a model replusive interaction with Vq = 20 o.u. and ro = 0.1 o.u. (dotted curve) 
corresponding to age — 0.01553 o.u. (224.3 Bohr). In the same figure, we also include 
the non-interacting (Vq = 0, a^c = 0) case (continuous curve), which naturally lies below 
the repulsive interaction (a^c > 0) curve. In Fig. 3, we plot ooQ{r) for an attractive 
interaction, viz., Vq = —100 o.u., ro = 0.0855 o.u (note from Fig. 1 that this corresponds 
to zero two-body bound state and age — —0.1176 o.u.) for A — 10. Since we cannot 
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go to large values of A due to numerical problems mentioned above, we keep A small 
and increase Vq to study the critical behaviour (see below) at a lower value of A. Both 
these curves have the general features same as those found in earlier calculations using 
K-harmonics approximation [13]. Fig. 3 shows a metastable region with a local minimum 
of u!o{r), which is preceded by a collapse region for smaller r. As A increases above a 
critical value (Acr), the metastable region disappears. This is seen in Fig. 4 for A = 16 
for the same Vq and tq. These features are the same as reported earlier [13]. However, 
in our case, since V{r) is finite for r — > 0, and the repulsive centrifugal term goes as ^, 
there is no real collapse. For very small r, u!o{r) becomes strongly repulsive even for an 
attractive two body interaction. This is represented by the dotted lines in Figs. 3 and 
4. Note that the dotted and continuous parts together constitute the entire calculated 
ujQ{r) curve. The small r (repulsive) part is plotted with a different (dotted) curve to 
emphasize that the remaining part (continuous portion) of uJo{r) has the same behaviour 
as obtained with attractive contact interaction in ref. [13]. Only the dotted part differs 
remarkably from the corresponding part in ref. [13]. In reality for A > A^r, there is a 
very narrow and deep well at a small value of r; hence all the particles will be trapped 
within this well. As the particles come within a small region, corresponding to a small 
value of r, the density of the condensate increases, and due to increased three and higher 
body collisions, molecule formation takes place with the disappearance of the BEG. The 
deep and narrow well in u!o{r) near the origin, for an attractive two-body interaction with 
A > Acr, can support a lowlying, highly localized bound state, which describes the for- 
mation of molecules. Although this is the lowest lying state in the corresponding a;o(r), 
it docs not represent the ground state of the condensate, which has already "collapsed". 
This gives a realistic scenario of what happens as A increases above Acr for attractive 
interactions. For an attractive (5-f unction interaction, the lack of a rigorous solution fails 
to give a reahstic picture and one talks of a "collapse of the condensate" in a quahtative 
fashion. 

We next calculate first three excited states for different number of par- 
ticles (A) in the condensate. These are shown in Fig. 5. Values of E^^ for n = 1, 2, 3 have 
been represented by diamonds, pluses and squares respectively. The excitation energy 
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increases slowly with A. They agree fairly well with the K-harmonic approximation [13]. 
In Table 2, wc present numerical values and notice that the excitation energies increase 
gradually with A. 

In Fig. 6, we plot the ground state wave function, Co{'>^), as a function 
of the global length r for various values of A. It is seen that as the particle number 
increases, the peak of Co('^) shifts towards larger values of r. This is understandable, 
since for large A, the total repulsion of all the pairs increases as A^ and particles are 
pushed outwards, by the A-dependent hypercentrifugal repulsion in eq. (28). 

Finally we calculate and plot the ground state energy per particle 
{Eo/A — ^fiw) as a function of Aagc for selected values of A (10, 20 and 30) for a repulsive 
interaction in Fig. 7. Corresponding curves are from the bottom upwards respectively. 
We compare these with the corresponding values calculated from the GP equation. This 
curve is the top most in Fig. 7. One notices that our results approach the GP result as A 
increases for a fixed Ausc, as expected. We also note that our energies are below those of 
the GP equation, indicating once again a better result from the variational point of view. 
Fig. 7 agrees qualitatively with a similar figure of ref. [19], where exact diagonalization 
of the Hamiltonian was performed for one and two dimensional condensates respectively. 

IV. Conclusions 

In this communication, we have investigated the T — properties of a 
Bosc-Eiiistciu condensate (BEG), consisting of A atoms (bosons) trapped by an external 
field and interacting via realistic two-body interactions. An ab initio treatment of the 
Schrodinger equation involves 3(^4 — 1) degrees of freedom for the relative motion. Use of 
traditional hyperspherical harmonics expansion (HHE) method is impossible for A > 3, 
due to tremendous and mounting complexity of the method as particle number increases 
beyond three. We circumvent this difficulty by exploiting the subset of potential har- 
monics (PM) basis, instead of the full set of hyperspherical harmonics (HH) basis. The 
PH basis is obtained as the subset of HH needed for expanding the two-body potential 
for the interacting pair. The choice of PH basis corresponds to inclusion of two-body 
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correlations and disregard of all higher-body correlations in the condensate. On the other 
hand, two-body correlations are very important in BEC and cannot be disregarded as 
in mean field theories or the GP equation. This assumption is exactly appropriate for 
the BEC, since for practical realization of BEC, the density of atoms must be kept so 
low that there are practically no three and higher body coUisions. Existence of the latter 
type of collisions would facilitate formation of molecules and consequent depletion of the 
condensate. As a consequence of this assumption, only four active degrees of freedom of 
the condensate (instead of a total of 3A — 3 degrees of freedom for the relative motion 
of the A particle system) are physically important - these are constituted by the global 
length (hyperradius, r) and the three active angle variables of the PH. In effect one freezes 
the remaining {3A — 7) angle variables of PH. This leads to a tremendous simplification 
of the actual numerical calculation. Since we make Faddeev like decomposition of the 
full wave function, an appropriate symmetrization of the wave function under exchange 
of the interacting pair guarantees full symmetrization. Moreover, the potential matrix 
elements involve integrals over only three angle variables, leading to an immense reduc- 
tion in the complexity of the numerical procedure for A. Since there are no theoretical 
restrictions on A, this opens the possibility of an approximate but very reliable, ab initio 
solution of the large but finite body condensate. However, a numerical difficulty arises 
due to the fact that the parameter a {— {3 A — 8)/2) of the Jacobi polynomials, P^'^{x), 
and its associated weight function, become very large as A increases. These cause nu- 
merical problems, for A > 40. We are at present attempting to remove this difficulty by 
appropriate numerical procedure. In the present report, we restrict ourselves to ^4 < 35, 
for which reliable calculations are possible. 

We have compared our results with earlier calculations for A — 3 [12], 
K-harmonic approximation [13], exact diagonalization of the Hamiltonian in one and two 
dimensions [19] as also with the predictions of the GP equation [8]. As a preliminary 
calculation we have taken two-body Gaussian interactions of varying range. Our results 
agree qualitatively with the previous ones, most of which use a contact interaction. This 
demonstrates the reliability and feasibility of our method. Thus a reliable ab initio calcu- 
lation for a large but finite number of atoms in a condensate, where individual particles 
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interact via realistic two-body interactions, appears feasible. Extension of our method to 
larger number of particles as also use of more realistic two-body interaction is underway. 
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Appendix 



Hyperspherical variables and hyperspherical 

harmonics 

Al : Hyperspherical variables 

The relative motion of the A — (N + 1) particle system is described in 
terms of N Jacobi coordinates defined by eq. (8) and having 3A^ degrees of freedom. An 
equivalent set of hyperspherical variables is constituted by the hyperradius (r) defined by 

eq.(ll), spherical polar angles of Ci, C2 , Ca^ and {N — 1) hyperangles {02, ^3, ■■■■4'n} 

giving the length of the Jacobi vectors Ci, C2 Xn, through 



Cn ^ r cos(f)N 
Cn-1 — f sin(f)N cos(f)N-i 
Cn-2 — f sin(f)N sm0jv-i cos(j)N-2 



(01 = 0) 

Eq. (39) automatically satisfies eq. (11). 



(39) 



C2 — sin(f)N sirKpN-i-'-sirKps cos(f)2 
Ci — r sin(t)N sirKpN-i-'-sirKps sin(l)2 



A2. Grand orbital operator 

The general grand orbital operator, L^(Ojv) of eq. (15) is defined 
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through [14] 

Lm - ^ + [3(^-2W4 + 2M0,-W.)]4- + £^ + %^ 

= 4(1 - zf)i^ + 6[2 - i{l + z,)]i: + + (40) 

(^ = 2, N) 



where Zi — cos2(f)i, uji represents the set of two polar angles of and 0j's are given by 
eq. (39). Note that L\{Q,i) — P{oui) and L%{^In) = L^{^n) appear in eq. (15). 



A3. Hyperspherical harmonics 

An eigenfunction of L^(Qjv) is called hyperspherical harmonics (HH) 
and is given (without angular momentum coupling) by [34] 

TV 

Yic]{nM) = Yi.mM n Yi^^^{u;,Y^^P'^f^-\<j>,) (41) 

J=2 



where 



{cos(t>j)^i (sm0j)^^-iPnT''''^^ (cos20,-) (i = 2,3,..., N) 



with 

Vj = z/j_i + 2nj + / j + I 

= £, + f-l 

^ ^ (43) 
= + 2nj + 

(j = 2,3,...,7V) 

In eq. (42) P^'^{x) is a Jacobi Polynomial. In eq. (41), [£] rcprcscnta the set of quantum 
numbers {(/i, mi), {I2, m2), (/at, mN),n2, n^, n^v} for a fixed value of grand orbital 
quantum number C — Cn- The quantum number d is defined through 

d = A_i + 2ni + /i (44) 
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with Ci — li. Hence 

N 

C = Cn = Ii + $](2n,- + /,) (45) 

The HH of eq. (41) forms the uncoupled basis. For systems with a good orbital angular 
momemtum L — Ix -\- I2 -\- ... -\- In^ one has to couple the individual orbital angular 
momenta - then the projection quantum numbers mi, m2, itlm are replaced by the 
[N — 1) intermediately coupled angular momenta and the projection M of L. 

The potential harmonics (PH) given by eq. (17) corresponds to Ijq — I, 
h — h — h — ■■■ — In-1 — 0, such that L — In — I, M — mjv = m and grand orbital JC 
= jCn — 2ir + 1 with n2 — — ... — un-i — and K — un. Substitution of these in 
eqs. (41) - (43) gives the PH of eq. (17). 
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Table 1. Calculated ground state energy per particle (in o.u.) of the condensate con- 
taining A particles for various K^ax values, showing convergence trend as Kmax increases 
{Vq = 3.1985 X 10^ o.u. and ro = 0.005 o.u.) 



A 


2 


4 


6 


8 


10 


12 


3 


1.50041 


1.50031 


1.50026 


1.50023 


1.50021 


1.50019 


5 


1.50123 


1.50117 


1.50112 


1.50108 


1.50104 


1.50101 


10 


1.50350 


1.50348 


1.50346 


1.50344 


1.50342 


1.50340 


15 


1.50453 


1.50451 


1.50450 


1.50449 


1.50449 




20 


1.50539 


1.50538 


1.50537 


1.50536 


1.50536 




25 


1.50618 


1.50617 


1.50617 


1.50616 


1.50616 




30 


1.50693 


1.50692 


1.50692 


1.50691 






35 


1.50761 


1.50763 


1.50763 


1.50766 







Table 2. Calculated excitation energies (in o.u.) of the first three excited states for dif- 
ferent numbers {A) of ^'^Rh atoms (parameters as in Table. 1). 



A 






3rd 


3 


2.00116 


4.00283 


6.00494 


5 


2.00130 


4.00428 


6.00962 


10 


2.00231 


4.00705 


6.01268 


15 


2.00355 


4.0130 


6.03147 


20 


2.00471 


4.01647 


6.04604 


25 


2.00671 


4.03225 


6.12762 


30 


2.03276 


4.08846 


6.17127 


35 


2.08319 


4.10225 


6.27118 



30 



-1000 -800 -600 -400 -200 200 400 

Vq {o.u.) 

Fig. 1 - Plot of calculated function of Vq for ro = 0.0855 o.u.. The dotted line 

corresponds to the Born approximation (a^f-*). 
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2 4 6 8 10 12 14 16 18 20 22 



r 

2 - Lowest eigen potential for ^4 = 20 as a function of r. Continuous curve is for 
(no two-body interaction) and the dotted curve is for a repulsive interaction 
0.01553 O.U.). 
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Fig. 3 - Plot of uJo{r) as a function of r (dotted and continuous curves together) for A 
(subcritical number) for a model attractive two-body interaction (Vq = —100 o.u., 
0.0855 O.U.), which corresponds to age — —0.1176 o.u. 



33 



200 



-200 
-400 [: 
Uo{r) -600 
-800 h 
-1000 -7 
-1200 h 
-1400 







10 



15 
r 



20 



25 



30 



Fig. 4 - Plot of uJo{r) as a function of r (dotted and continuous curves together) for A 
(critical number) for the same attractive two-body potential as in Fig. 3. 
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Fig. 5 - Three low- lying excitation frequencies for various values of particle number {A), 
corresponding to the JILA experiment with ^^Rb atoms ( a^c = 100 Bohr, oscillator fre- 
quency = 200 Hz). Energies are in oscillator units. 
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Fig. 6 - Plot of ground state wave function as a function of hyperradius (r) for various 
indicated value of A, in the chosen trap. 
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Fig. 7 - Plot of ground state energy per particle {Eq/A — 3/2huj) as a function of Aag 
for a repulsive interaction for indicated values of A and the GP results. 
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